# Compute the power spectra of detrended and non-detrended, and compare
RR = 1./(max(rv_dict['jd'])-rv_dict['jd'][0])
df = RR/2
f_N = 1./(2*np.median(np.diff(rv_dict['jd']))) # Nyquist frequency
f = np.arange(0, f_N, df)
GLS_object_trended = LombScargle(time, RV_offset_corrected, normalization = 'psd') # Pre-plan the periodogram
pls_trended = GLS_object_trended.power(f[1:]) # Exclude 0 frequency to avoid divide by zero
GLS_object_detrended = LombScargle(time, RV_offset_corrected-GF, normalization = 'psd') # Pre-plan the periodogram
pls_detrended = GLS_object_detrended.power(f[1:]) # Exclude 0 frequency to avoid divide by zero
#================ Same thing but now remove planet b. Can we find planet c??
orbit_b = kep_model.keplerian['0'].rv(kep_model.t)
GLS_object_trended_bremoved = LombScargle(time, RV_offset_corrected-orbit_b, normalization = 'psd') # Pre-plan the periodogram
pls_trended_bremoved = GLS_object_trended_bremoved.power(f[1:]) # Exclude 0 frequency to avoid divide by zero
GLS_object_detrended_bremoved = LombScargle(time, RV_offset_corrected-GF-orbit_b, normalization = 'psd') # Pre-plan the periodogram
pls_detrended_bremoved = GLS_object_detrended_bremoved.power(f[1:]) # Exclude 0 frequency to avoid divide by zero
#=============
fig, axes = plt.subplots(1, 2, sharey=True, sharex=True, figsize = (9, 7), dpi = 250.)
axes[0].loglog(f[1:], normalize(RV_offset_corrected, f, pls_trended), '-', color = "red", label = "No detrending", alpha = 0.5)
axes[0].loglog(f[1:], normalize(RV_offset_corrected-GF, f, pls_detrended), 'b-', label = "GF detrending", alpha = 0.5)
axes[1].loglog(f[1:], normalize(RV_offset_corrected-orbit_b, f, pls_trended_bremoved), '-', color = "red", label = None, alpha = 0.5)
axes[1].loglog(f[1:], normalize(RV_offset_corrected-GF-orbit_b, f, pls_detrended_bremoved), 'b-', label = None, alpha = 0.5)
axes[0].text(0.5*RR, 1.1*10**6, "$\mathcal{2R}$")
axes[0].fill_betweenx(np.linspace(10**1, 10**6, 1000), 2*RR, x2 = 0, color = 'k', alpha = 0.2)
axes[0].vlines([1./5574.2], 10**1, 10**6.1, color = 'r', alpha = 0.8, ls = '--', label = "55 Cnc d", zorder = 0)
axes[0].vlines(f_mag, 10**1, 10**6.1, color = 'g', alpha = 0.8, ls = '--', label = "Mag. Cycle", zorder = 0)
axes[1].text(0.5*RR, 1.1*10**6, "$\mathcal{2R}$")
axes[1].fill_betweenx(np.linspace(10**1, 10**6, 1000), 2*RR, x2 = 0, color = 'k', alpha = 0.2)
axes[1].vlines([1./5574.2], 10**1, 10**6.1, color = 'r', alpha = 0.8, ls = '--', label = None, zorder = 0)
axes[1].vlines(f_mag , 10**1, 10**6.1, color = 'g', alpha = 0.8, ls = '--', label = None, zorder = 0)
axes[1].scatter(lunar, 10**5.15, s = 15.0, marker = 'o', color = 'black', label = 'Moon')
axes[1].scatter(2./len(time), 10**5.3, s = 15.0, marker = 'o', color = 'cyan', label = '$f = 2/N$')
axes[0].scatter(planet_b, 10**6.2, s = 15.0, marker = 'o', color = 'orange', label = '55 Cnc b')
axes[1].scatter(planet_c, 10**4.61, s = 15.0, marker = 'o', color = 'magenta', label = '55 Cnc c', zorder = 0)
axes[0].legend(loc = 'lower left', fontsize = 'small')
axes[1].legend(loc = 'lower left', fontsize = 'small')
axes[0].annotate('$\mathcal{3R/4}$', xy=(np.mean([f_mag, (1./5574.2)]), 10**6.23), xytext=(np.mean([f_mag, (1./5574.2)]), 10**6.28), xycoords='data',
ha='center', va='bottom',
bbox=dict(boxstyle='square', fc='white', color='white'),
arrowprops=dict(arrowstyle='-[, widthB=0.58, lengthB=0.8', color='k'))
axes[1].annotate('$\mathcal{3R/4}$', xy=(np.mean([f_mag, (1./5574.2)]), 10**6.23), xytext=(np.mean([f_mag, (1./5574.2)]), 10**6.28), xycoords='data',
ha='center', va='bottom',
bbox=dict(boxstyle='square', fc='white', color='white'),
arrowprops=dict(arrowstyle='-[, widthB=0.58, lengthB=0.8', color='k'))
axes[0].set_ylim(bottom = 10**1); axes[1].set_ylim(bottom = 10**1)
axes[0].set_xlim(left = f[1], right = 10**(-1)); axes[1].set_xlim(left = f[1], right = 10**(-1))
axes[0].set_xlabel("Frequency [cycles/day]")
axes[1].set_xlabel("Frequency [cycles/day]")
axes[0].set_ylabel("Normalized Power")
axes[0].set_title("No Planets Removed")
axes[1].set_title("55 Cnc b removed")
plt.subplots_adjust(wspace = 0.05)
plt.show()
